
package org.zhongyin.Regression;
import java.sql.Timestamp;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import org.zhongyin.Dao.LocalDBDaoImpl;

import Jama.Matrix;

/**
 * 
 * @author ZhongYin Zhang
 *
 */

public class Regression {
    
    private static final double TUNING_CONSTANT = 4.685d;
    
    private Matrix coefficient;
    
    private double leastSquares;
    
    private double standardDeviation;
    
    private Matrix weight;
    
    private static final boolean DEBUG = false;
    
    //D = sqrt(eps(class(double))) from STATROBUSTFIT_Revised.m under matlab
    private static final double D = 1.490116119384766e-08;
    
    private static final int iterLimit = 50;
    
    public Regression(){
            coefficient = null;
            
            leastSquares = 0d;
            
            standardDeviation = 0d;
            
            weight = null;
            
    }
    
    
    
    /**
     * For Linear Regression
     * b = (XT*X)-1 * XT * Y
     * @param array_X       X
     * @param array_Y       Y
     * @return      coefficient (b)
     */
    public Matrix getCoefficient(double[][] array_X, double[] array_Y){
            
            if(array_X == null || array_Y == null){
                    throw new NullPointerException();
            }
            Matrix X = new Matrix(array_X);
            Matrix Y = new Matrix(array_Y, array_Y.length);
            
            Matrix X_Transpose = X.transpose();
            // XT*X
            Matrix X_Transpose_Times_X = X_Transpose.times(X);
            //(XT*X)-1 Matrix inverse
            Matrix square_inverse = X_Transpose_Times_X.inverse();
            //(XT*X)-1 * XT
            Matrix inverseTimesTranspose = square_inverse.times(X_Transpose);
            
            Matrix result = inverseTimesTranspose.times(Y);
            
            
            coefficient = result;
            
            return result;
            
    }
    
    
    // for linear regression
    /**
     * @param array_Y       Y
     * @param array_X       X
     * @param array_b       coefficient
     * @return      leastsquares (e2)
     */
    public double getLeastSquares(double[] array_Y, double[][] array_X, double[][] array_b){
            if(array_Y == null || array_X == null || array_b == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            Matrix b = new Matrix(array_b);
            Matrix Xb = X.times(b);
            
            
            double[][] array_Xb = Xb.getArray();
            
            double[][] e = new double[array_Y.length][1];
            
            
            double result = 0d;
            
            
            for(int i = 0; i < e.length; i++){
                    
                    e[i][0] = array_Y[i] - array_Xb[i][0];
                    result += e[i][0] * e[i][0];
            }
            
            
            
            leastSquares = result;
            
            return result;
            
            
    }
    
    
    /**
     * ri = (yi - ybari) / std
     * @param array_Y       Y
     * @param array_X       X
     * @param array_b       coefficient
     * @param sDeviation 
     * @return      Standardized residuals
     */
    public Matrix getStdResidual(double[] array_Y, double[][] array_X, double[][] array_b, double sDeviation){
            if(array_Y == null || array_X == null || array_b == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            Matrix b = new Matrix(array_b);
            Matrix Xb = X.times(b);
            
            double[][] y_bar = Xb.getArray();
            
            double[] r = new double[array_Y.length];
            
            for(int i = 0; i < r.length; i++){
                    r[i] = (array_Y[i] - y_bar[i][0])/sDeviation;
            }
            
            Matrix stdResidual = new Matrix(r, array_Y.length);
            
            return stdResidual;
    }
    
    
    
    /**
     * 
     * Compute the adjusted residuals
     * @param array_Y Y
     * @param array_X X
     * @param array_b coefficient
     * @param array_h leverage
     * @return
     */
    public double[] getAdjustedResiduals(double[] array_Y, double[][] array_X, double[][] array_b, Matrix array_h){
            if(array_Y == null || array_X == null || array_b == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            Matrix b = new Matrix(array_b);
            Matrix Xb = X.times(b);
            
            double[][] h = array_h.getArray();
            
            double[][] y_bar = Xb.getArray();
            
            double[] r = new double[array_Y.length];
            
            for(int i = 0; i < r.length; i++){
//                  System.out.print(array_Y[i] - y_bar[i][0]+ ", ");
                    r[i] = (array_Y[i] - y_bar[i][0])/Math.sqrt(1-h[i][i]);
//                  r[i] = 1/Math.sqrt(1-h[i][i]);
            }
            
//          System.out.println();
            
            return r;
    
    }
    
    /**
     * s is the robust variance given by MAD/0.6745 where MAD 
     * is the median absolute deviation of the residuals
     * 
     * @param array_Radj  Radj
     * @param array_X X
     * @return
     */
    public double getS(double[] array_Radj, double[][] array_X){
            if(array_Radj == null){
                    return 0;
            }
            
            
             double[] array_Radj_Clone = array_Radj.clone();
            
             double median = 0d;
             
         int length = array_X[0].length; 
         
//       System.out.println ("length : " + length);
         for(int i = 0; i < array_Radj.length; i++){
             array_Radj_Clone[i] = Math.abs(array_Radj[i]);
         }
            
             Arrays.sort(array_Radj_Clone);
             
             
             
             array_Radj_Clone = Arrays.copyOfRange(array_Radj_Clone, length-1, array_Radj.length);
             
//           System.out.println("Sorted Radj: " + Arrays.toString(array_Radj_Clone));
             
             if(array_Radj.length % 2 == 0){
                     
                     
                     median = (array_Radj_Clone[array_Radj_Clone.length/2] + array_Radj_Clone[array_Radj_Clone.length/2 - 1])/2;
             }
             else{
                     median = array_Radj_Clone[array_Radj_Clone.length/2];
             }
             
//           median = array_Radj[array_Radj.length/2];
             
//           System.out.println("Media :" + median);
             double MAD = median;
             
             double S = MAD / 0.6745;
             
             return S;
            
            
    }
    
    /**
     * 
     * @param stdResidual ri
     * @return      W[][] with Wii(ri)
     */
    public Matrix getOriginalWeight(Matrix stdResidual){
            if(stdResidual == null){
                    throw new NullPointerException();
            }
            
            double[][] r = stdResidual.getArray();
            
            double[][] W = new double[r.length][r.length];
            
            for(int i = 0; i < r.length; i++){
                    
                    if(Math.abs(r[i][0]) <= TUNING_CONSTANT){
                            W[i][i] = (1 - (r[i][0] * r[i][0]) / (TUNING_CONSTANT * TUNING_CONSTANT))*(1 - (r[i][0] * r[i][0]) / (TUNING_CONSTANT * TUNING_CONSTANT));
                    }
                    else{
                            W[i][i] = 0d;
                    }
            }
            
            
            
            Matrix Weight = new Matrix(W);
            
            weight = Weight;
            
            return Weight;
            
            
    }
    
    public Matrix getRobustWeight(Matrix U){
            if(U == null){
                    throw new NullPointerException();
            }
            
            double[][] u = U.getArray();
            
            double[][] W = new double[u.length][u.length];
            
            for(int i = 0; i < W.length; i++){
                    
                    if(Math.abs(u[i][0]) <= 1){
                            W[i][i] = (1 - (u[i][0] * u[i][0]))*(1 - (u[i][0] * u[i][0]));
                    }
                    else{
                            W[i][i] = 0d;
                    }
            }
            
            
            
            Matrix Weight = new Matrix(W);
            
            weight = Weight;
            
            return Weight;
            
            
    }
    
    
    /**
     * 
     * The quantities hi are the leverages
     * , the diagonal elements of X(Xt*X)-1 * Xt
     * @param array_X X
     * @return array_h
     */
    public Matrix getLeverage(double[][] array_X){
            if(array_X == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            
            // h = X(XTX)-1XT
            Matrix result = X.times((X.transpose().times(X)).inverse()).times(X.transpose());
            
            return result;
    }
    
    
    /**
     * @param leastSquares e2
     * @return e
     */
    public double getStandardDeviation(double leastSquares){
            
            double result = Math.sqrt(leastSquares);
            
            standardDeviation = result;
            
            return result;
            
    }
    
    
    
    /**
     * For Weighted Linear Regression
     * b = (XT*WX)-1 * XT * W *Y
     * @param array_X       X
     * @param array_Y       Y
     * @return      coefficient (b)
     */
    public Matrix getWeightedCoefficient(double[] array_Y, double[][] array_X, Matrix weight){
            if(array_X == null || array_Y == null || weight == null){
                    throw new NullPointerException();
            }
            Matrix X = new Matrix(array_X);
            Matrix Y = new Matrix(array_Y, array_Y.length);
            
            Matrix X_Transpose = X.transpose();
            
            //XT*W
            Matrix X_Transpose_Times_W = X_Transpose.times(weight);
            // XT*W*X
            Matrix X_Transpose_Times_W_Times_X = X_Transpose_Times_W.times(X);
            //(XT*WX)-1 Matrix inverse
            Matrix square_inverse = X_Transpose_Times_W_Times_X.inverse();
            //(XT*WX)-1 * XT
            Matrix inverseTimesTranspose = square_inverse.times(X_Transpose);
            //(XT*WX)-1 * XTW
            Matrix timesWeight = inverseTimesTranspose.times(weight);
            //(XT*WX)-1 * XTWY
            Matrix result = timesWeight.times(Y);
            
            
            coefficient = result;
            
            return result;
    }
    
    
    
    
    /**
     * For Weighted Linear Regression
     * @param array_Y       Y
     * @param array_X       X
     * @param array_b       coefficient
     * @return      leastsquares (e2)
     */
    public double getWeightedLeastSquares(double[] array_Y, double[][] array_X, double[][] array_b, Matrix weight){
            if(array_Y == null || array_X == null || array_b == null || weight == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            Matrix b = new Matrix(array_b);
            Matrix Xb = X.times(b);
            
            
            double[][] array_Xb = Xb.getArray();
            
            double[][] e = new double[array_Y.length][1];
            
            
            double result = 0d;
            
            
            for(int i = 0; i < e.length; i++){
                    
                    e[i][0] = array_Y[i] - array_Xb[i][0];
                    //Wii*(Yi-Xib);
                    result += weight.getArray()[i][i] * (e[i][0] * e[i][0]);
            }
            
            
            
            leastSquares = result;
            
            return result;
            
            
    }
    
    
    /**
     * wfun in Matlab
     * @param array_Radj
     * @param S
     * @return
     */
    public Matrix getU(double[] array_Radj, double S){
            if(array_Radj == null){
                    return null;
            }
            
            double[] u = new double[array_Radj.length];
            for(int i = 0; i < array_Radj.length; i++){
                    u[i] = array_Radj[i]/(TUNING_CONSTANT*S);
            }
            
            Matrix U = new Matrix(u, u.length);
            
            return U;
            
            
    }
    /**
     * Iteratively get the weight by Irls
     * @param array_Y Y
     * @param array_X X
     * @param weight array of weight
     * @return get more accurate weight
     */
    public Matrix getIterWeight(double[] array_Y, double[][] array_X, Matrix weight, int type){
            
//          Regression regression = new Regression();
            
            int length = weight.getArray().length;
            
            int iter = 0;
            
//          //Generate an identity matrix
//          weight = Matrix.identity(length, length);
            
            double[] weight_array = new double[length];
            
            for(int i = 0 ; i < length; i++){
                    weight_array[i] = weight.get(i, i);
            }
            

            // First coefficient (b)
            Matrix array_b = getWeightedCoefficient(array_Y, array_X, weight);

            // First h
            Matrix leverage = getLeverage(array_X);

            // First adjusted residuals
            double[] array_Radj = getAdjustedResiduals(array_Y, array_X, array_b.getArray(), leverage);

            // First S sigma in matlab
            double S = getS(array_Radj, array_X); 

            // First U
            Matrix U = getU(array_Radj, S);

            
            
            if(DEBUG){
              System.out.println(iter + " weight: " + Arrays.toString(weight_array));
              System.out.println(iter + " b: " + Arrays.deepToString(array_b.getArray()));
              System.out.println(iter + " h: " + Arrays.deepToString(leverage.getArray()));
              System.out.println(iter + " Radj: " + Arrays.toString(array_Radj));
              System.out.println(iter + " S: " + S);
              System.out.println(iter + " U: " + Arrays.deepToString(U.getArray()));
              System.out.println();
            }
            
            
            
//          
            
            
            
            while(any(array_b.getArray()) && iter < iterLimit){
                    
                    iter++;
                    
                    weight = getRobustWeight(U);
                    
                    for(int i = 0 ; i < length; i++){
                            weight_array[i] = weight.get(i, i);
                    }
                    

                    
                    array_b = getWeightedCoefficient(array_Y, array_X, weight);

                    // First leverage
                    leverage = getLeverage(array_X);

                    // First adjusted residuals
                    array_Radj = getAdjustedResiduals(array_Y, array_X, array_b.getArray(), leverage);

                    // First S
                    S = getS(array_Radj, array_X); 

                    // First U
                    U = getU(array_Radj, S);

                    if(DEBUG){
                      System.out.println(iter + " weight: " + Arrays.toString(weight_array));
                      System.out.println(iter + " b: " + Arrays.deepToString(array_b.getArray()));
                      System.out.println(iter + " h: " + Arrays.deepToString(leverage.getArray()));
                      System.out.println(iter + " Radj: " + Arrays.toString(array_Radj));
                      System.out.println(iter + " S: " + S);
                      System.out.println(iter + " U: " + Arrays.deepToString(U.getArray()));
                      System.out.println();
                    }
                    
            }
            
            
            double[][] diagonal_weight = new double[weight_array.length][weight_array.length];
            
            for(int i = 0 ; i < weight_array.length; i++){
                    diagonal_weight[i][i] = weight_array[i];
            }
            Matrix result = new Matrix(diagonal_weight);
            
            
            
            
            double RMSE =  getRobustMSE(array_Radj, array_X, S, leverage);
            
            LocalDBDaoImpl l = new LocalDBDaoImpl();
            l.getEntityManager().getTransaction().begin();
    		
            l.addRMSE(type, RMSE);
    		
    		l.getEntityManager().flush();
    		l.getEntityManager().getTransaction().commit();
//    		l.getEntityManager().close();
            
            System.out.println("RMSE : " + RMSE);
            return result;
    }
    
    /**
     * 
     * True if any element of a vector is
     * logical 1 (TRUE).
     * @param array
     * @return
     */
    public boolean any(double[][] array){
            if(array == null){
                    return false;
            }
            
            for(int i = 0; i < array.length; i++){
                    if(Math.abs(array[i][0])  <= D * Math.max(Math.abs(array[i][0]), 0)){
                            return false;
                    }
            }
            
            return true;
    }
    
    
    /**
     * get forgetting factor based on lumda from 0.970 to  1.000
     * @param weight
     * @return a list of forgetting factor
     */
    public List<double[][]> getWeightPrimePrime(Matrix weight){
            if(weight == null){
                    throw new NullPointerException();
            }
            
            double k = 0.971d;
            
            List<double[][]> list = new ArrayList<double[][]>();
            
            int length = weight.getArray().length;
            
            
            while( k <= 0.971){
                    Matrix lumda = Matrix.identity(length, length);
                    double[][] array_lumda = lumda.getArray();
                    
                    System.out.print("w: ");
                    for(int i = 0; i < length; i++){
                            array_lumda[i][i] = lumda.get(i, i) * Math.pow(k, length-1 - i) * weight.get(i, 0);
                            System.out.print(array_lumda[i][i]+ ", ");
                    }
                    

                    list.add(array_lumda);
                    k += 0.001;
            }
            
            return list;
    }
    
    public double getSecondStageWeightedLeastSquare(double[] array_Y, double[][] array_X, List<double[][]> weightPrimePrime){
            if(array_Y == null || array_X == null || weightPrimePrime == null){
                    throw new NullPointerException();
            }
            
            double minE = Double.MAX_VALUE;
            double minLeastSquare = Double.MAX_VALUE;
            
            int counter = 0;
            for(double[][] array_weight: weightPrimePrime){
                    
                    if(counter == 29){
                            System.out.println(counter + "weight: " + Arrays.deepToString(array_weight));
                    }
                    
//                  System.out.println(counter + "weight: " + Arrays.deepToString(array_weight));
                    
                    Matrix weight = new Matrix(array_weight);
                    Matrix coefficient = getWeightedCoefficient(array_Y, array_X, weight);
                    
                    double[][] array_b = coefficient.getArray();
                    
//                    System.out.println(counter+ "b:" + Arrays.deepToString(array_b));
                    
                    double weightedLeastSquare = getWeightedLeastSquares(array_Y, array_X, array_b, weight);
                    
                    double e = Math.sqrt(weightedLeastSquare/array_Y.length);
                    
//                    System.out.println(counter+ "e2: " + weightedLeastSquare);
//                    
//                    System.out.println(counter+ "e: " + e);
                    
                    minE = Math.min(minE, e);
                    
                    minLeastSquare = Math.min(minLeastSquare, weightedLeastSquare);
                    
                    counter++;
                    
            }
            
            return minE;
    }
    
    
    public List<Matrix> getFFCofficient(double[] array_Y, double[][] array_X, List<double[][]> weightPrimePrime){
            if(array_Y == null || array_X == null || weightPrimePrime == null){
                    throw new NullPointerException();
            }
            
            double lumbda = 0.971d;
            
            
            
            List<Matrix> result = new ArrayList<Matrix>(); 
            for(double[][] array_weight: weightPrimePrime){
                    
                    
                    

                    Matrix weight = new Matrix(array_weight);
                    Matrix coefficient = getWeightedCoefficient(array_Y, array_X, weight);
                    
                    System.out.println();
                    
                    System.out.println(lumbda + " b: " + Arrays.deepToString(coefficient.getArray()));
                    
                    result.add(coefficient);
                    
                    
                    lumbda += 0.001;
                    
            }
            
            return result;
            
            
    }
	
    /**
     * Compute the robust MSE, translate function statrobustsigma(wfun,r,p,s,t,h) in matlab
     * @param array_Radj Radj
     * @param array_X
     * @param S sigma
     * @param leverage h
     */
    public double getRobustMSE(double[] array_Radj, double[][] array_X, double S, Matrix leverage){
    	
    	// r in matlab
    	Matrix Radj = new Matrix(array_Radj, array_Radj.length);
    	//xrank
    	int p = array_X[0].length;
    	//n= length(r)
    	int length = array_Radj.length;
    	// st = s*t
    	double st = S*TUNING_CONSTANT;
    	
    	
    	double[] array_st = new double[length];
    	for(int i = 0; i < array_st.length;i++){
    		array_st[i] = st;
    	}
    	Matrix m_st = new Matrix(array_st, array_st.length);
//    	Radj =Radj.inverse();
    	// u = r./st;
    	Matrix m_u = Radj.arrayRightDivide(m_st);
    	
    	
    	//feval(wfun, u)
    	Matrix m_w = getRobustWeight(m_u);
    	double[] array_w = new double[length];
    	for(int i = 0; i < array_w.length; i++){
    		array_w[i] = m_w.getArray()[i][i];
    	}
    	
    	m_w = new Matrix(array_w,array_w.length);
    	
    	//phi = u.*feval(wfun, u)
    	Matrix phi = m_u.arrayTimes(m_w);
    	
    	double delta= 0.0001;
    	double[] array_delata = new double[length];
    	
    	for(int i = 0; i < array_delata.length; i++){
    		array_delata[i] = delta;
    	}
    	Matrix m_delta = new Matrix(array_delata, array_delata.length);
    	//u1 = u-delta;
    	Matrix m_u1 = m_u.minus(m_delta);
    	
        m_w = getRobustWeight(m_u1);
    	 array_w = new double[length];
    	for(int i = 0; i < array_w.length; i++){
    		array_w[i] = m_w.getArray()[i][i];
    	}
    	
    	m_w = new Matrix(array_w,array_w.length);
    	
    	//phi0 = u1 .* feval(wfun,u1);
    	Matrix phi0 = m_u1.arrayTimes(m_w);
    	//u1 = u + delta;
    	m_u1 = m_u.plus(m_delta);
    	
    	m_w = getRobustWeight(m_u1);
	   	array_w = new double[length];
	   	for(int i = 0; i < array_w.length; i++){
	   		array_w[i] = m_w.getArray()[i][i];
	   	}
	   	
	   	m_w = new Matrix(array_w,array_w.length);
    	
    	
    	//phi1 = u1 .* feval(wfun,u1);
    	Matrix phi1 = m_u1.arrayTimes(m_w);
    	
    	//dphi = (phi1 - phi0) ./ (2*delta);
    	Matrix dphi = (phi1.minus(phi0)).arrayRightDivide(m_delta.times(2));
    	
    	
    	double[] array_dphi = dphi.getColumnPackedCopy();
    	
    	double sum = 0d;
    	for(int i =0; i < array_dphi.length;i++){
    		sum+=array_dphi[i];
    	}
    	
    	double mean = sum/array_delata.length;
    	
    	//m1 = mean(dphi)
    	double m1 = mean;
    	
    	double[] array_h = new double[length];
    	double[] array_1 = new double[length];
    	for(int i = 0; i < array_h.length; i++){
    		array_h[i] = leverage.getArray()[i][i];
    		array_1[i] =1;
    	}
    	
    	Matrix h = new Matrix(array_h, array_h.length);
    	Matrix constant1 = new Matrix(array_1, array_1.length);
    	
    	Matrix constant1Minush = constant1.minus(h);
    	
    	double[] array_phi = phi.getColumnPackedCopy();
    	for(int i = 0; i < array_phi.length;i++){
    		array_phi[i] = array_phi[i]*array_phi[i];
    	}
    	
    	Matrix phiSquare = new Matrix(array_phi, array_phi.length);
    	
    	Matrix temp = constant1Minush.arrayTimes(phiSquare);
    	
    	sum = 0;
    	
    	for(int i = 0 ; i < length; i++){
    		sum+=temp.get(i, 0);
    	}
    	
    	//m2 = sum((1-h).*phi.^2)/(n-p);
    	double m2 = sum/(length - p);
    	//K = 1 + (p/n) * (1-m1) / m1;
    	double K = 1 + (p/length) * (1-m1) / m1;
    	
    	//s = K*sqrt(m2) * st /(m1);
    	double result = K * Math.sqrt(m2) * st / m1;
    	
    	return result;
    	
    }
    
    
    
    

	
	
	
}
